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1.  Introduction. 

Many  applications  in  science  and  enginering  involve  mixed  systems  of  differential  and 
algebraic  equations  (DAE’s).  For  some  examples  see,  for  instance,  the  monograph  [BCPS9]. 
It  is  hardly  surprising  that  such  systems  share  many  properties  with  ordinary  differential 
equations  (ODE’s).  In  fact,  recent  existence  theories  ([Rh91],  [RR91a],  [RR91b])  have 
shown  that,  in  general,  a  DAE  can  be  reduced  locally  to  an  (explicit)  ODE  on  some 
submanifold  of  the  space  of  unknown  variables. 

However,  despite  the  strong  analogy  between  DAE’s  and  ODE’s,  important  differences 
exist.  For  instance,  from  the  fact  that  DAE’s  are  reducible  to  ODE’s  only  on  some  sub¬ 
manifold  of  the  solution  space  it  follows  that  solutions  of  a  DAE  can  pass  only  through 
points  on  such  a  submanifold;  that  is,  its  initial  values  must  satisfy  certain  compatibility 
conditions.  Beyond  this,  solutions  of  DAE’s  may  exhibit  features  that  solutions  of  explicit 
ODE’s  cannot  possess.  For  instance,  the  simple  problem 

xj  +  X2  =  0,  X2  =  1.  z(0)  =  (1,  -1) 

has  the  unique  solution  x(i)  =  ((1  —  1 11/2,  t  —  1)  which  cannot  be  continued  beyond  t  =  1 
despite  the  fact  that  ,r(l)  =  (0,0)  and  limf—1-  x(t)  =  a’(l)  exist.  This  situation  would  be 
impossible  for  solutions  of  explicit  ODE’s. 

'The  work  was  supported  in  part  by  ONR-grant  N-00014-90-J-1025,  NSF-grant  CCR-8907654.  and 
AFOSR-grant  900094. 

2Department  of  Mathematics  and  Statistics.  University  of  Pittsburgh,  Pittsburgh,  PA  15260 
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In  the  electrical  engineering  literature  such  type  of  points  have  been  called  impasse 
points  (see  e.g.  [C69]  or  [CDS9]  where  also  other  references  are  given).  Although  they 
have  no  analog  in  connection  with  explicit  ODE’s  they  are  closely  related  to  the  “singular 
points”  of  implicit  ODE’s.  In  [Ra89]  the  most  often  encountered  type  of  such  singularities 
for  implicit  ODE’s  was  analyzed,  called  there  standard  singular  points.  In  a  recent  paper 
[RR92a]  it  was  shown  that  the  geometric  reduction  theory  for  DAE’s  presented  in  [RR91b] 
allows  for  a  generalization  of  the  results  in  [Ra89]  to  so-called  standard  impasse  points  of 
quasilinear  DAE’s. 

The  aim  of  this  article  will  be  to  show  that  the  theory  of  [RR92]  leads  naturally  to 
the  development  of  a  computational  procedure  for  the  explicit  computation  of  standard 
impasse  points  of  quasilinear  DAE’s.  For  this  we  outline,  in  Sections  2  and  3,  briefly  and 
without  proof,  some  of  the  relevant  results  for  singular  ODE’s  and  DAE’s  from  the  cited 
earlier  papers.  Then  Section  4  presents  the  details  of  the  computational  algorithm  and 
finally  in  Section  5  we  give  some  numerical  examples  which  show  the  effectivity  of  the 
process. 

2.  Singular  Points  of  ODE’s. 

Definition  2.1.  Consider  a  quasilinear  problem 

(2-1)  B(y)y  =  H(y),  y(0)  =  y0 

where  B  :  V  — *  £(Rn)  and  H  :  V  — >  Rn  are  C 1  on  some  open  set  V  C  Rn.  A  point  y  €  V 
is  a  regular  point  of  (2.1)  if  rank  B(y)  =  n  and  a  singular  point  if  rank  B(y)  <  n  but  y  is 
a  limit  point  of  regular  points  of  (2.1 ). 

Clearly,  for  a  regular  point  y0  €  V  the  initial  value  problem  (2.1)  has  a  unique  solution 
in  a  neighborhood  of  yo-  But  already  simple  examples  show  that  the  behavior  of  the 
solutions  of  (2.1)  in  a  neighborhood  of  a  singular  point  yo  may  vary  strongly  with  the 
type  of  singularity  encountered  there.  A  partial  classification  of  singularities  which  will  be 
sufficient  for  our  purposes  is  given  next: 
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Definition  2.2.  (i)  A  singular  point  y  6  V  of  the  ODE  in  (2.1)  r-singular  if 

(2.2)  dimkerB(y)  =  r. 

(ii)  A  1 -singular  point  is  called  basic  if 

(2.3)  H{y)i  rgeB(y). 

(iii)  A  basic  1 -singular  point  is  a  standard  singular  point  if 

(2.4)  DB(y)(u,u)  £  rge  B(y),  Vu  6  kerB(y)  \  {0}. 

We  summarize  here  briefly  the  theory  developed  in  [Ra89]  for  the  case  of  standard 
singular  points  (see  also  [RR92]).  With  a  standard  singular  point  y0  as  starting  point  the 
initial  value  problem  (2.1)  cannot  have  a  C1  solution  y  :  J  — ►  V  on  an  open  interval  J 
containing  the  origin.  In  fact,  this  would  require  that  £(yo)y(0)  =  H(y0)  which  contradicts 
(2.3).  Thus  at  a  standard  singular  point  we  may  expect  at  best  “one-sided”  solutions  in 
the  following  sense: 

Definition  2.3.  With  a  standard  singular  point  yo  €  V  as  starting  point  a  solution  of 
the  initial  value  problem  (2.1)  is  any  continuous  function  y  \  J  —>  T>  defined  on  an  interval 
J  =  [0,  T)  or  J  —  (— T,  0]  for  some  T  >  0  which  is  of  class  C1  on  J°  =  J  \  {0}  and  satisfies 
2/(0)  =  y0  and  B(y(t))y(t)  =  H(y(t))  fort  <E  J°. 

With 


(2.5)  a(y)(u,  v)  =  (vTH(y))  ( vTDB(y)(u , «)) 

the  two  conditions  (2.3)  and  (2.4)  for  a  standard  singular  point  y  are  equivalent  to 

(2.6)  a(y)(u,v)  ^  0,  Vu  €  kerB(y)  \  {0},  v  €  ker B(y)T  n  {0}. 
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Since  (2.5)  is  a  continuous,  quadratic  form  in  u  and  in  v  its  value  must  be  either  positive 
or  negative  for  all  pairs  of  nonzero  vectors  u  £  ker  B(y),  v  £  kevB(y)T  if  only  this  holds 
for  one  such  pair  of  vectors. 

The  principal  existence  result  for  solutions  near  standard  singular  points  can  now  be 
phrased  as  follows  (see  [Ra89,  Theorem  5.1])  : 

Theorem  2.1.  Let  y0  £  V  be  a  standard  singular  point  of  the  ODE  in  (2.1).  Then  the 
initial  value  problem  (2.1)  has  exactly  two  solutions  which  a le  both  defined  on  J  =  [0,  T) 
or  on  J  =  (— T,  0]  for  some  T  >  0  depending  upon  a(y0)(u,v )  >  0  or  ct(y0)(u,v)  <  0, 
respectively,  for  some  pair  of  nonzero  vectors  (u,  v)  £  kerB(y)  x  ker  B(y)r.  Moreover, 
||y(f)||  tends  to  infinity  as  t  £  J  \  {0}  tends  to  zero. 

Theorem  2.1  implies  that  a  solution  of  (2.1)  starting  at  some  regular  point  can  reach 
a  standard  singular  point  yo  at  some  later  time  only  if  the  form  (2.5)  is  negative  at  y°. 
Standard  singular  points  y0  with  positive  o(y0)  obviously  can  never  be  reached  in  increasing 
time.  Thus,  in  view  of  the  theorem,  the  following  notation  is  appropriate: 

Definition  2.4.  A  standard  singular  point  yo  of  (2.1)  is  accessible  or  inaccessible  if 
Q(yo)(w,u)  <  0  or  a(y0)(u,  v)  >  0,  respectively,  for  some  pair  of  nonzero  vectors 

{u,v)  £  ker  .B(y)  x  ker  l?(y)T. 

The  theorem  asserts  that  accessible  standard  singular  points  are  reached  in  finite  time 
by  trajectories  emanating  elsewhere  in  V.  Since  these  trajectories  cannot  be  continuously 
extended  beyond  these  points,  they  represent  “catastrophes”  for  the  solutions  of  (2.1)  and 
standard  ODE-solvers  fail  near  such  points.  It  can  also  be  shown  (loc.  cit.)  that  no  small 
perturbation  of  the  initial  condition  (and/or  of  B  or  H )  can  affect  the  eventual  encounter 
of  such  points. 

Standard  singular  points  are  analogous  to  limit  points  of  parametrized  nonlinear  equa¬ 
tions 

F(z,  A)  =  0. 

Suppose,  indeed,  that  F  :  V  — ►  R"  is  of  class  C1  on  some  open  set  V  C  Rn  x  R1  and 
that  z  :  J  — >  R"  is  a  C1  mapping  on  an  open  interval  J  such  that  (r(A),A)  £  V  and 
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F(z(\),  A)  =  0  for  A  £  J.  Then 


(2.7) 


DzF(z(  A),  A) 

0 


0\  (z'( A)\  /  D\F(z(\),  A) 

JUM  i 


where  primes  indicate  differentiation  with  respect  to  A,  and  hence  (2.7)  is  an  ODE  of 
the  quasi-linear  form  (2.1).  The  singular  points  of  (2.7)  are  exactly  those  (z,A)  6  T>  for 
which  rank  DzF(z ,  A)  <  n;  that  is,  which  are  foldpoints  of  F  with  respect  to  A  (see  e.g. 
[Rh86]).  The  simplest  foldpoints  are  the  limit  points  for  which  dimker  DzF(z,  A)  =  1  and 
D\F(z,  A)  ^  rge  DzF(z,  A).  Obviously,  these  two  properties  correspond  to  the  conditions 
(2.2)  and  (2.3)  and  it  is  readily  checked  that  (2.4)  holds  exactly  for  the  simple  limit  points 
of  F  (with  respect  to  A),  (see  e.g.[Rh86]). 

Foldpoints  of  a  parametrized  nonlinear  system  are  typically  computed  by  solving  a 
suitably  augmented  form  of  the  system.  It  is  natural  to  consider  the  same  approach  for 
the  computation  of  singular  points  of  ODE’s. 

Let  y  =  y(t)  be  a  C^-solution  of  (2.1)  and  suppose  that  a  C1  function  r  :  R1  — >  R1  with 
strictly  positive  derivative  is  used  to  define  a  transformation  t  —  r(s )  of  the  independent 
variable  t.  Then  rj(s )  =  y(r(s ))  satisfies 


BMts  =  TsHM- 


By  Theorem  2.1  the  derivative  dy/dt  becomes  infinite  when  the  solution  approaches  a 
standard  singular  point  y* .  This  suggest  that  r  should  be  chosen  such  that  dr/ds  tends  to 
zero  as  we  approach  y*  but  drj/ds  remains  bounded.  For  instance,  we  may  wish  to  specify 
r  implicitly  by  using  a  normalization  cT(drj/ds )  =  1  with  a  suitable  vector  c  €  R". 

This  normalization  may  be  obtained  by  means  of  an  augmented  system  of  the  form 


(2.8) 


fB{y)  /v(t/)\  /0\ 

UT  o  JWy))~  W 


where  c  €  K"  is  chosen  such  that  at  a  point  y  €  T>  under  consideration  the  matrix  of  (2.8) 
is  nonsingular.  Certainly,  such  an  augmentation  can  be  found  if  and  only  if  y  is  either  a 
regular  point  of  (2.1)  or  a  basic  1-singular  point. 
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Hence  the  matrix  of  (2.8)  remains  nonsingular  for  all  y  in  some  open  neighborhood 
U  C  T>  of  y  whence  for  fixed  y  €  U  the  solution  (u(y),7(y))  €  Rn+1  of  (2.8)  is  unique. 
Obviously  we  have  v(y)  ^  0  for  all  y  £  U  and  7 (y)  ^  0  for  all  regular  points  y  €  £/  of  (2.1). 
Moreover,  because  of  (2.3)  we  see  that  7(y*)  =  0  at  any  basic  1-singular  point  y*  €  ZA 

For  any  regular  point  yo  G  ZY  the  initial  value  problem  (2.1)  has  a  unique,  local  C1- 
solution  y  :  [0,T)  — >  ZY  for  some  T  >  0.  Suppose  that  7(yo)  >  0.  Then  y(y(t))  >  0  for 
0  <  t  <  T  and  the  initial  value  problem 

(2.9)  §  =  7(y(<)),  *  €  [0,  T),  t(0)  =  0, 

has  a  unique,  monotonically  increasing  solution  r  :  [0,  cr)  — >  R1  with  0  <  a  <  T.  Hence, 
as  desired,  r  defines  a  transformation  of  the  independent  variable  t  of  (2.1).  As  before,  we 
set 


(2.10)  t?(s)  =  y(r(s)),  s€[0,<r),  77(0)  =  y0- 


Together  with  (2.9)  the  chain  rule  provides  that 


(2.11) 


=  7(77(5))  ^(r(s)),  s  €  [0,  cr) 


whence  by  (2.1)  and  (2.8)  it  follows  that 


B{y{t))  -  v(y(<))]  =  0,  0  <  t  <  T\ 


and  therefore,  by  (2.10)  and  (2.11),  that 

(2.12)  ^(s)  =  v(rj(s)),  rj(0)  =  y0,  -s  €  [0, a). 
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Proposition  2.1.  Suppose  that  the  solution  y  =  y(t)  of  (2.1)  tends  to  some  standard 
singular  point  y*  €  U  and  hence  has  been  extended  to  a  maximal  interval  [0,  T*)  such 
that  y(t )  =  y*.  (Necessarily,  y*  is  a  basic  1-singular  point  since  the  matrix  of  the 

augmentation  (2.8)  is  invertible  at  y*  by  hypothesis).  Then  (i)  there  exists  cr*  <  oc  such 
that  the  solution  t  =  r(s)  of  (2.9)  is  defined  for  s  €  [0,<7*)  and  that  lims—^*  r(s)  =  T* . 
(ii)  The  solution  of  (2.12)  is  defined  and  of  class  C 1  on  [0,  cr*  +  e)  for  some  e  >  0.  (Hi)  If, 
in  addition  y*  is  a  standard  singular  point,  then  7(77(5))  changes  sign  as  s  crosses  <7*. 

Proof:  Let  J  denote  the  set  of  all  a  such  that  there  exists  a  C1  -solution  of  (2.9)  for 
s  €  [0,  a)  satisfying  0  <  t(s)  <  T*  on  that  interval.  Clearly,  J  is  not  empty  and  hence 
a*  =  sup{  a  :  cr  6  J)  is  well  defined.  Thus  there  exists  a  ^-solution  of  (2.9)  for  5  6  [0,  cr*) 
satisfying  0  <  t(s )  <  T*  on  this  interval.  In  order  to  show  that  lims_CT.  r(s)  =  T*  note 
that  r  remains  monotonically  increasing  on  [0,  cr*)  and  hence  that  lim,_(T.  r(s)  =  r*  <  T* 
exists.  Suppose  that  r*  <  T*.  Then,  by  the  continuity  of  7,  the  compactness  of  y([0,  r*]), 
and  the  fact  that  7 (y(t))  >  0  for  t  €  [0,  T*),  it  follows  that  there  exists  a  positive  constant 
7o  for  which  7 (y(t))  >  70  in  [0,X*).  This  implies  that  0*  <  00,  for  otherwise  r(s)  >  705 
for  5  €  [0,oo)  and  hence  lims_oo  r(s)  =  00  in  contradiction  with  0  <  t($ )  <  r*  <  00  for 
s  E  [0,  cr*)  =  [0, 00).  But  now,  by  setting  r(a*)  =  r*,  we  can  define  a  continuous  extension 
of  r  to  some  interval  [0,cr*  +  e]  with  sufficiently  small  e  >  0  such  that  0  <  r(s)  <  T*. 
Thus,  assuming  r*  <  T*  we  obtain  a  contradiction  with  the  maximality  of  cr*.  This  shows 
that  r*  =  T*. 

The  above  arguments  show  only  that  a*  <  00,  but  we  can  show  that  cr*  <  00.  In  fact, 
assume,  to  the  contrary,  that  <7*  =  00  so  that  77  is  defined  in  [0, 00).  Since  by  construction 
cTv(y(s))  =  1  it  follows  from  (2.12)  that  cT (dr) / ds)(s)  =  1  whence  cTr/(s)  =  cTy0  +  s.  This 
implies  that  lim3_00 cTrj(s)  =  oc  and  hence  also  that  lima_oo  ||77(s)j|  =  00.  But  then  we 
arrive  at  a  contradiction  since  ((77(^7) ff  =  |jy(7'(s))||  and  ||y(f)||  is  bounded  on  the  compact 
interval  [0,T*].  This  proves  (i). 

It  follows  from  (i)  that  (2.12)  has  a  unique  solution  77  for  5  6  [0, <7*)  and  rj(<r*)  =  y* 
defines  a  continuous  extension  of  77.  Therefore,  the  solution  of  (2.12)  can  be  extended  to  a 
larger  interval  [0,cr*  4-  e)  with  some  c  >  0  as  claimed  in  (ii). 


As  noted  earlier,  we  have  7(77(3))  >  0  for  0  <  s  <  a*  and  y(rj(cr*))  =  7 (y*)  =  0.  Thus 
a*  is  the  first  zero  of  7(77(3))  in  [0,<7*  +  e).  If  now  y*  is  a  standard  singular  point,  7 (77(3)) 
must  change  sign  as  3  crosses  a*.  For  this  note  that  by  differentiation  of 

B(r}(s))v(r](s))  =  ■j(t](s))H(tj(s)),  s  £  [0,ct*], 
together  with  7 (y*)  =  0  we  obtain  for  all  h  £  R"  that 

DB(y*)(h,  v(y*))  +  B(y*)Dv(y')h  =  (D^{y*)h)H{y*). 


For  h  =  v(y*)  ^  0  and  any  nonzero  vector  w  £  ker  B{y*)T  it  follows  from  (2.5)  and  (2.6) 


^  wTDB(y*)(v(y*),v(y*))  a(y*)(v{y*),w) 

Dl(y  My)  = - - =  WWrW '  # 

Since  ^(7 orj)/ds)(a*)  =  D^{y*)v{y*),  this  proves  that  indeed  (7077 )(s))  must  change  sign 
as  s  crosses  a*.  □ 

Altogether,  therefore,  by  solving  (2.12)  and  monitoring  the  first  sign  change  of  7(77(3)) 
we  can  calculate  a*  and  hence  y*  =  77(0-*).  The  value  of  T*  is  then  given  by 


(2.13)  T*  =  f  7(77(3))*, 

Jo 

which  follows  directly  from  (2.9)  and  lim,—^*  77(3)  =  T*. 

The  augmentation  procedure  described  here  is  designed  to  work  in  the  neighborhood 
of  a  standard  singular  point.  But  in  practice,  also  higher  order  singularities  y*  £  V  are 
encountered  where  the  matrix  of  the  simple  augmented  system  (2.8)  becomes  singular.  In 
order  to  avoid  difficulties  near  such  points,  we  may  work  with  an  overdetermined  augmented 
system  of  the  form 


/B{y)  —H(y)  E\ 

V  cT  0  0  J 


/  V(y)  \ 

I  w{y)T  j 

\Z(y)TJ 


(2.14) 
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Here,  for  given  q,  1  <  q  <  n,  the  matrices  E  and  C  have  dimension  n  x  (q  —  1)  and  r?  x  q, 
respectively,  and.  correspondingly,  in  the  solution.  V(y),  w(y).  and  Z(y)  are  blocks  of  size 
n  x  q,  q  x  1,  and  q  x  (q  —  1),  respectively.  As  before,  the  matrices  E  and  C  are  chosen  such 
that  at  some  “current”  point  yc  €  T>  the  matrix  of  (2.14)  is  nonsingular  and  hence  remains 
nonsingular  for  all  y  in  some  open  neighborhood  Uc  C  T>  of  yc.  Thus,  for  each  y  6  Uc 
the  solution  of  (2.14)  is  unique. Clearly,  for  sufficiently  large  q  this  can  be  accomplished 
irrespective  of  dim  ker  B(yc )  and  even  in  the  case  when  H(yc)  €  rge  B{yc ):  (except  when 
H(yc)  —  0;  but  see  Remark  2.1  further  below). 

We  summarize  some  basic  properties  of  the  augmentation  (2.14): 

Proposition  2.2.  For  given  q  >  1  and  y  €UC  the  solution  of  (2.14)  satisfies 

(2.15)  dim  ker  B(y )  =  dim  ker  (w(y),  Z(y)), 

and  rank  Z(y)  —  q  —  1  implies  that  rank  (B(y), —H(y))  =  n.  The  converse  holds  if 
H(y)  $  rge  B(y)  (and  hence  y  is  a  basic  1 -singular  point). 

Proof:  Generally,  for  y  €.  Uc  we  have 

(2.16)  B(y)V(y)  =  (ff(  y),  -E)(w(y).  Z(y)f, 

as  well  as  rank  V(y)  =  q  (since  CTV(y)  =  7?)  and  rank  (7f(y),  —E)  =  q  (since  the  matrix 
of  (2.14)  is  invertible)  which  together  imply  the  first  assertion.  Indeed,  since  both  V(y) 
and  ( H(y ),  —  E)  are  n  x  q,  q  <  n,  and  have  maximum  rank  q ,  we  infer  from  (2.16)  that 

dim  ker  B(y)V(y)  =  dim  ker  (w(y),  Z(y))T  =  dim  ker  (u'(y),  Z(y)^ 

(recall  that  (w(y),  Z(y))  is  q  x  q).  Moreover,  using  again  the  fact  that  ker  V(y)  =  {0}.  we 
see  that 

(2.17)  dim  ker  B(y)V(y)  <  dim  ker  B(y). 
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Thus,  to  complete  the  proof  of  (2.15),  it  suffices  to  show  that  the  converse  inequality  of 
(2.17)  holds.  This  follows  at  once  if  we  can  show  that  ker  B(y)  C  rge  V(y).  For  this  let 
u  £  ker  B(y ),  so  that 


(2.18) 


'B(y)  —H(y)  E 

CT  0  0 


On  xl 

CTU 


u 

Oi  Xl 

‘  0(9  —  1)  Xl 

On  the  other  hand,  CTV(y)  =  Iq  implies  that  CtV(T)Ctu  =  CTu  and  hence,  by  (2.16), 
that 


(2.19) 


(B(y)  - H(y )  E\ 

\  cT  o  o ) 


/  V(y)CTu  \ 

I  w(y)TCTu  I 

V  Z(y)TCTu  / 


But  the  systems  (2.18)  and  (2.19)  have  the  same  right  hand  side,  and  hence,  because  the 
matrix  is  invertible,  the  solutions  are  identical  whence,  in  particular.  ?/  =  V(y)CTu  and 
therefore  u  £  rge  V(y). 

If  rank  Z(y)  =  q  —  1  then  rge  EZT  —  rge  E.  Therefore,  rge  E  C  rge  ( B(y),—H(y )) 
and  hc-m  rge  (B(y),  -H{y),E)  =  rge  (B(y),  -H(y))  which  for  rank  ~H(y))  <  u 

contradicts  the  nonsingularity  of  the  matrix  of  (2.14). 

To  prove  that,  conversely,  rank  (B(y),  —H(y))  =  n  and  H(y)  £  rge  B(y)  imply 
rank  Z(y)  —  q  —  1  suppose  that  H(y)  ^  rge  B(y)  and  rank  Z(y)  <  q  —  2.  so  that 
dim  ker  Z(y)T  >  2.  Let  a  =  1,2  be  two  linearly  independent  vectors  in  ker  Z(y)T. 
By  (2.16)  we  have 

B(y)V(y)u0  =  ( w{y)TuQ)H{y ),  o  =  1.2 

and  hence  w{y)T uQ  =  0,  a  =  1,2,  since  H(y )  £  rge  B(y).  Thus  uQ  £  ker  (zr(y).  Z(y))1 . 
o  =  1,  2.  By  (2.15)  the  linear  independence  of  the  two  vectors  implies  that  dim  ker  B(y)  > 
2.  which  in  turn  implies  that  rank  ( B(y ),  —H(y))  <  n.  □ 

For  y  £UC  and  any  nonzero  vector  a(y)  £  ker  Z(y)T  we  have 


(2.20o)  B(y)v(y)  =  ~/(y)H(y),  c(y)Tv(tj)  =  a{y)T  a{y). 

(2.20  b)  v(y)  =  V(y)a(y).  c(y)  =  Co(y).  ->(y)  =  ir{y)Ta(y). 
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This  has  the  general  form  of  (2.8),  and  as  before,  j(y)  =  0  implies  that  B{y)  is  singular. 
But,  the  converse  is  not  necessarily  true.  Moreover,  the  vector  n(ij)  must  depend  smoothly 
on  y.  This  is  easily  guaranteed  as  long  as  rank  Z(y )  —  q  —  1  but  is  not  generally  feasible 
unless  we  drop  the  assumption  that  a(y )  ^  0.  Let  Zl(y)  denote  the  (q  —  1)  x  (q  —  1) 
submatrix  obtained  from  Z(y)  by  deleting  the  z-th  column.  Then  the  vector 

(2.21)  a(y)  =  (ai(y),...  ,aq(y))7’,  a{(y)  =  (-l)‘det  Z'(y),  i  =  1, . . .  ,y, 

obviously  depends  smoothly  on  y  and  satisfies  Z(y)Ta(y)  =  0,  (see  e.g.  [S73],  Appendix 
II).  Moreover,  we  have  a(y)  ^  0  exactly  if  rank  Z(y)  =  q  —  1. 

For  this  choice  of  a(y)  the  following  result  holds: 

Proposition  2.3.  For  y  G  Uc  and  with  the  vectors  (2.21 )  consider  the  relations  (2.20). 
Then  7 (y)  =  0  exactly  if  B(ij)  is  singular. 

Proof:  If  7(y)  =  0  then  either  rank  Z(y)  <  q  —  1  in  which  case,  by  (2.15),  B(y)  is 
singular,  or  rank  Z(y)  =  q  —  1  whence  a(y)  ^  0  and  thus  also  v(y )  ^  0  which  together 
with  B(y)v(y)  =  0  implies  again  that  B(y)  is  singular.  Conversely,  suppose  that  B{y)  is 
singular.  In  the  case  rank  Z(y)  -  q  -  1,  we  have  again  v(y)  ^  0  and  Proposition  2.2 
ensures  that  H(y)  £  rge  B(y).  Thus  the  first  equation  (2.20a)  leads  to  a  contradiction 
unless  7(y)  =  0.  On  the  other  hand  for  rank  Z(y)  <  q  -  1  we  necessarily  have  7(y)  =  0 
because  o(y)  =  0  in  that  case.  □ 

By  Proposition  2.2  we  see  that  when  B(y )  is  singular  and  H(y )  G  rge  B(y)  then 
necessarily  rank  Z(y)  <  q  —  1.  Such  points  are  evidently  not  standard  impasse  points, 
and  hence  constitute  “higher”  singularities.  Our  choice  (2.21)  of  a(y)  evidently  transforms 
these  points  into  equilibrium  points  of  the  dynamic  system  (2.12)  with  v  as  in  (2.20b). 

Computationally  the  simplest  case  arises  with  q  =  2  where  (2.14)  has  the  form 


(B(y) 

-H(y)  e\ 

(My) 

My)  \ 

(° 

°\ 

T 

cl 

0 

0 

m(y) 

My)  = 

1 

0 

V  cj 

0 

0/ 

\  My) 

My) ' 

VO 

J 
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(2.22) 


Hence,  the  vector  (2.21)  becomes  here  a(y )  =  ( — ^2(2/),  -^i(y))  and  we  obtain 

v(y)  =  z2(y)vi(y )  -  2i(y)v2(y), 

(2-23)  7 (y)  =  z2(y)w1(y)  -  Zi(y)ut2(y), 

c(y)  =  z2{y)ci{y)  -  z\(y)c2(y). 

Proposition  2.3  implies  that  7 (y)  does  not  vanish  on  a  trajectory  terminating  at  a  singu¬ 
lar  point  y *  €  Uc.  Obviously  (2.22)  has  the  same  general  form  as  the  simple  augmentation 
(2.8)  and  the  computational  procedure  is  the  same  as  before,  namely,  we  form  and  solve 
the  explicit  equation  (2.12).  The  only  difference  is  that  the  normalization  condition  now 
involves  the  nonconstant  vector  c(y)  while  in  (2.8)  this  vector  was  constant. 

The  constancy  of  c  was  used  to  prove  that  for  a  basic  1-singular  point  there  exists  cr*  <  00 
such  that  the  solution  t  —  r(s)  of  (2.9)  is  defined  for  [0, <7*)  and  that  lim,—**  r(s)  =  T*. 
The  result  is  easily  extended  to  the  case  when  c  depends  on  y.  For  this  note  that  the 
proof  of  the  existence  of  a*  <  00  carries  over  verbatim  a-  !  irrespective  of  the  singularity 
encountered  at  y*;  all  that  is  needed  is  that  the  matrix  (2.22)  be  invervible  at  y*.  Now  if 
y*  is  a  basic  1-singular  point  and  <7*  =  00  we  conclude  from  c(r](s))Tv(rj(s))  =  ||a(y)||22 
that  c(y*)Tv(rj(s))  >  e  >  0  for  all  s  close  to  cr*.  Indeed,  rank  Z(y*)  =  q  —  1  by  Propo¬ 
sition  2.2  and  hence  a(y*)  ^  0.  In  other  words,  for  sufficiently  large  s  it  follows  that 
c(y*)T(dr]/ds)(s)  >  e  and  therefore  that  c(y*)Tr](s)  >  e(s  —  s0)  +  c(y*)Ty(s0)  for  s  >  s0 
and  so  large  enough.  As  before,  this  contradicts  the  boundednes  of  ||y(.s)||  for  s  >  0. 

If  y*  is  a  standard  singular  point,  then  7 (y(s))  changes  sign  as  s  crosses  a*.  Obviously, 
if  y*  is  not  a  basic  1-singular  point  then  we  can  no  longer  ascertain  that  a *  <  00. 

Remark  2.1  The  matrix  of  any  augmentation  (2.14)  and,  in  particular,  that  of  (2.8) 
or  (2.22),  will  be  singular  at  any  point  y  €  Uc  where  H(y)  =  0;  that  is,  at  any  stationary 
(equilibrium)  point  of  the  autonomous  ODE  (2.1).  This  reflects  the  fact  that  any  regular 
stationary  point  can  be  reached  only  in  infinite  time  and  hence  our  scale-transformation 
must  become  indetermined  along  such  trajectories.  Clearly,  the  natural  resolution  of  this 
difficulty  is  to  make  the  system  (2.1)  non- autonomous  by  adding  the  equation  t  =  1. 


3.  Singular  Points  of  DAE’s. 

In  this  section  we  turn  to  differential-algebraic  systems  of  the  following  form: 
Definition  3.1.  The  equation 

(3.1)  A(x)x  =  G(x), 

with  C2  mappings  A  :  T>  — >  £(R")  and  G  :  V  — *  Rn  on  some  open  set  T>  C  R"  is  a 
quasilinear  DAE  on  V  if 

(3.2)  Gr(x)  €  rge  A(x),  x  6  T>  =>  rank  A(x)  =  r 
and  if  the  mapping 

(3.3  )  (x,p)  €  V  x  Rn  i — >  A(x)p  -  G(x)  e  Rn, 

is  a  submersion. 

The  submersion  property  of  (3.3)  requires  that  for  every  (x,p)  €  Z>  x  R"  the  mapping 

(3.4)  (/i,Jt)€Rn  x  R”  i — *■  (DA(x)h)  p  +  A(x)k  -  DG(x)h  €  Rn, 
is  onto  (see  e.g.  [AMR88]).  As  a  consequence  the  set 

(3.5)  M  =  {(x,p)  eDxR”:  A(x)p  -  G(x)  =  0} 

is  a  closed  n-dimensional  C2 -submanifold  of  D  x  R". 

In  [RR91b]  a  geometric  procedure  was  developed  for  reducing  an  implicit  DAE  F(x,  x)  = 
0  to  a  system  of  ODE’s  on  a  manifold  locally  near  a  point  (x°,p°)  €  F_1(0).  A  simplified 
version  of  this  reduction  process  for  quasilinear  DAE’s  (3.1)  is  given  in  [RR92].  In  that 
case  the  reduction  is  local  only  in  the  first  variable  due  to  the  linearity  of  the  equation  in 
x.  We  summarize  briefly  this  process  for  (3.1). 
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Set 


(3.6)  W  =  {x  E  V  :  G(x )  €  rge  .A(x)}, 

so  that  ( x,p )  €  M  if  and  only  if  x  €  W  and  hence  W  =  tt(M)  where  7r  :  T>  x  Rn  — ►  R"  is 
the  projection  onto  the  first  factor.  Under  the  conditions  of  Definition  3.1  it  can  be  shown 
(see  [RR92],  Proposition  3.1)  that  W  is  an  r-dimensional  C2-submanifold  of  T>  and  that 
W  is  closed  in  V  if  the  set  {16P:  rank  A(x)  =  r}  is  closed  in  T>.. 

For  any  C1-solution  x  :  J  — ►  T>  of  (3.1)  on  an  open  interval  J  C  R1  we  must  have 
x(t)  €  W  for  t  €  J  and  thus  (x(t),x(t))  G  TW  for  all  t  6  J  where  we  view  here  the 
tangent  bundle  TW  as  a  subset  of  TRn  =  Rn  x  Rn.  Hence  (x(t),x(t))  €  M  implies  that 
( x(t),x(t ))  €  TW  fl  M  for  all  t  6  J.  The  desired  reduction  of  (3.1)  now  requires  a  (local) 
characterization  of  TW  fl  M. 

For  this  let  x°  6  W.  Then  there  exist  open  subsets  U  C  W  and  VcRr  and  a  C2-map 
<p  :  V  — *  Rn  which  is  a  diffeomorphism  of  V  onto  U.  In  other  words,  (f~1  is  a  chart  of  W 
at  x°. 

Evidently,  U  and  V  may  be  chosen  small  enough  such  that  there  is  a  linear  subspace 
Z°  C  Rn  which  complements  rge  A(<p(y))  for  all  y  €  V.  Let  P°  be  the  projection  of  Rn 
onto  rge  A(x(i)  along  Z°  and  L°  any  linear  isomorphism  from  rge  A(x°)  onto  Rr.  Then 
L°  o  P°  is  a  linear  isomorphism  of  rge  A(<p(y))  onto  Rr  for  all  y  €  V  and  it  follows  that 

f  x  =  <p(y),  p  =  D<p{y)q, 

{(x,p)  e  tw  r\  M,  x  eU)  < 

l  B(y)q  -  H(y)  =  0, 

where  we  have  set 

(3.7)  B(y)  =  L’pOAMy^D^y),  B(y)  =  L°P0GMy)). 

Evidently,  the  operators  B  and  H  map  into  £(Rr)  and  Rr,  respectively,  and  are  of  class 

C1. 
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If  the  interval  of  definition  J  C  R1  of  the  solution  is  restricted  to  ensure  that  x(  J)  C  U 
and  therefore  that  ( x(t),x(t ))  6  TW  fl  M  and  x(t)  €  U  for  all  f  6  7,  then  the  C1  function 

(3.8)  y  :  J  -*F,  y(t)  =  v>_1  o  x(<), 

is  a  C1 -solution  of  the  equation 

(3-9)  B{y)y  =  H(y), 

called  the  reduction  of  (3.1)  near  x°.  Conversely,  for  any  C^-solution  y  :  J  — *  V  of  (3.9) 
the  function  x(t)  —  ip  o  y(t)  is  a  C^-solution  of  (3.1). 

Evidently,  if  B(y)  is  invertible  in  a  neighborhood  of  y°  =  <,j-1(x0)  then  (3.9)  satisfies  the 
conditions  of  Definition  2.1  and  the  augmentation  procedure  of  Section  2  can  be  applied. 
This  is  the  case  when  (3.1)  has  index  1  in  the  sense  of  the  following  definition: 

Definition  3.2.  The  quasilinear  system  (3.1 )  is  a  nonsingular  DAE  of  index  1  if 

(3.10)  {x  e  W,  G(x)  €  rge  A(x)|rilv  }  =»  rank  A(x)\tzw  —  rank  A(x)(=  r). 

From  (3.10)  it  follows  that  B(y°)  €  <j£(Rr,  rge  A(x0))  and  D<p{y°)  €  GL(Rr,  TxoW ) 
and  this  provides  the  basis  of  the  following  existence  and  uniqueness  theorem  for  (3.1): 

Theorem  3.1.  Let  (3.1)  be  a  nonsingular  DAE  of  index  1.  Then,  for  any  x°  €  W\  = 
7 t(TW  fl  M)  C  W  there  exists  a  unique  C1  solution  x  :  J  — ►  V  on  some  open  interval  J 
containing  the  origin,  of  the  initial  value  problem 

(3.11)  A(x)x  =  G(x),  x(0)  =  x°. 

Moreover,  no  C1  solution  of  (3.11 )  exists  for  x°  £  Wj. 

Definition  3.2  does  not  rule  out  the  existence  of  points  x°  €  W  where  rank  A(x°)|ji0  w  < 
r.  Such  points  do  not  belong  to  the  set  Wi  =  n(TW DM)  and  hence  no  C1  solution  to  the 
corresponding  initial  value  problem  (3.11)  exists.  Nevertheless,  in  analogy  to  Definition 
2.3  ‘one-sided’  solutions  may  well  occur  at  such  points: 
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Definition  3.3.  A  solution  of  the  initial  value  problem  (3.11)  at  a  point  x°  ^  tt(TW  D  M) 
is  any  continuous  function  x  :  J  —*  V  defined  on  an  interval  of  the  form  J  =  [0,  T)  or 
J  =  (— T,  0]  for  some  T  >  0  which  is  of  class  C 1  on  J°  =  J  \  {0}  and  satisfies  x(t*)  =  x * 
and  A(x{t))x(t)  =  G(x(t))  for  t  €  J°. 


In  [RR92]  a  precise  definition  of  accessible  and  inaccessible  impasse  points  of  nonsingular 
DAE’s  of  index  1  (and  higher)  is  given  where  the  existence  of  one-sided  solutions  can  be 
guaranteed.  We  shall  not  repeat  this  theory  here  but  summarize  the  main  result  in  the 
form  of  the  following  theorem:  (see  [RR92],  Lemma  5.1  and  Theorem  5.1): 

Theorem  3.2.  Let  (3.1)  be  a  nonsingular  DAE  with  index  1.  The  point  x°  €  W  is  an 
accessible  or  inaccessible  impasse  point  of  (3.1 )  if  and  only  if  £0  =  cp-1  (x°)  is  an  accessible 
or  inaccessible  standard  singular  point,  respectively,  of  the  reduction  (3.9)  of  (3.1 )  locally 
near  x° .  Then,  the  initial  value  problem  (3.11)  has  exactly  two  solutions  in  the  sense  of 
Definition  3.3,  both  defined  either  on  J  —  [0,  T)  or  J  —  (-T,  0]  for  some  T  >  0.  Moreover, 
||i(t)||  tends  to  infinity  as  t  tends  to  zero. 


Because  C^-solutions  of  a  (not  necessarily  nonsingular)  DAE  (3.1)  lie  in  W\  —  it(W  D 
TM),  their  closure  relative  to  the  open  set  V  C  Kr  must  lie  in  W  when  W  is  closed  in 
V.  This  is  the  case  in  many  practical  applications  but,  mathematically,  it  is  not  the  only 
possibility.  When  W  is  not  closed  in  V  it  becomes  possible  for  points  x°  ^  W  to  be 
reached  in  finite  time  by  C1 -trajectories  that  cannot  be  continuously  extended  beyond 
x°.  Likewise,  there  may  be  points  x°  6  W  corresponding  to  higher  singularities  of  the 
reduction;  that  is,  with  dimker5(y°)  >  1,  at  which  Cl  trajectories  stop. 

As  an  illustration  we  consider  the  first  and  third  example  of  [CD89]  which  have  the  form 


(1  0  ox 
0  1  0 
\0  0  0  J 


X  =  Gj(x),  X  e  R3,  j  -  1,2 


with 


Gi(x)  = 


x3 

*3 


G2(x)  = 


Xi  A  x3(x2  +  x\)J 


16 


Hence,  in  both  examples  the  conditions  of  Definition  3.1  are  satisfied  with  r  =  2.  Moreover, 
we  have 

W  =  {x  €  R3  :  a?i  +  x3(r2  +  x2)  =  0}, 
and,  evidently,  the  mapping 


<^:R2  — »  R3,  ip(y)  = 


(  V2  Vi  +  3 y\ 


Dtp(y)  = 


-1 

\  0 


0 

-1 


is  here  a  global  diffeomorphism  from  R2  onto  W .  Therefore,  with  the  linear  isomorphism 


L°P° 


0  0\ 

1  0/’ 


the  reductions  of  the  two  problems  have  the  form 


/  V2  yi+  3 y\ 

\-l  0 


V  ~  Hj(y),  j  =  1,2 


where 


Hi(y)  = 


H2{y)  = 


y\ 


~V2 


Obviously,  in  both  cases,  the  singular  points  form  the  one-dimensional  manifold 


K  =  {y  €  R2  :  t/i  +  3 y\  =  0} 


and  we  have  dim  ker  B(y*)  =  1  and  rge  B(y *)  =  span  (yj,—  1)T  for  y*  G  K.  Moreover, 
for  j  =  1  we  see  that  Hi(y*)  £  rge  B(y*)  for  all  y*  €  K  with  y£  ^  —  1  while  for  j  =  2  we 
have  Hi(y*)  €  rge  B(y*)  for  all  y*  6  A\  Thus,  in  the  second  case,  none  of  the  points  of 
K  is  a  standard  impasse  point.  On  the  other  hand,  for  j  =  1,  a  simple  calculation  with 
u  =  (0,1)T  €  ker  B(y*)  and  v  =  (l,y2)T  €  ker  B(y*)T  shows  that 

a(y*)(u,v)  =  -6(y2)2(y2  +  1), 
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whence,  all  points  y*  €  I\  with  y|  ^  0,-1  are  here  standard  impasse  points  which  are 
accessible  for  y\  >  —  1,  y\  ^  0  and  inaccessible  for  y\  <  —  1. 

Thus,  for  j  =  1  the  points  y *  €  K,  y%  =  0,-1  and  for  j  =  1  all  points  y*  €  K  are 
higher  singularities.  In  both  cases,  y*  =  0  is  also  a  stationary  point.  The  differences 
between  these  higher  singularities  and  the  standard  impasse  points  in  the  first  example 
were  observed  in  [CD 89]  but  explained  differently. 


4.  Computation  of  Singular  Points  of  DAE’s. 

The  reduction  process  for  quasilinear  DAE’s  sketched  in  the  previous  section  and  the 
resulting  theory  of  impasse  points  for  nonsingular  DAE’s  of  index  1  suggests  that  we 
may  compute  such  points  by  applying  the  augmentation  approach  of  Section  2  to  the 
reduced  system  (3.9).  In  this  section  we  show  that  this  does  indeed  lead  to  an  efficient 
computational  algorithm. 

For  simplicity  we  develop  the  method  only  for  DAE’s  of  the  following  form  occurring 
frequently  in  applications: 


(4.1) 


^(x)’ 

0 


x  = 


Gi(x)' 

G2{x)i 


Here  A%  :  T>  — ►  £(Rr),  G\  :  V  —*  Rr,  G2  ’■  T>  — >  Rp  are  C2-maps  on  some  open  set  V  C  R", 
with  n  =  r  +  p,  p  >  0,  and  (4.1)  is  assumed  to  be  a  quasilinear  DAE  in  the  sense  of 
Definition  3.1. 


The  submersion  condition  for  the  mapping  (3.3)  requires  that  rank  DG2(x)  =  p  for 
x  6  T>  and  hence  that  the  set 


(4.2)  N  =  {xeV:  G2(x)  =  0} 

is  an  r-dimensional  C2 -submanifold  of  S. 

The  manifolds  (3.5)  and  (3.6)  are  here  given  by 

M  =  {(x,p)  €  V  x  Rn  :  x  €  N,  A,{x)p  =  Gj(x)} 
W  =  {x  €  Ar  :  rank  Ai(x)  =  r). 


(4.3) 

(4.4) 
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and  this  allows  for  a  considerable  simplification  of  the  reduction  process.  As  in  [RhS6]  we 
introduce  at  any  “current”  point  xc  £  N  a  tangential  local  coordinate  system.  For  this  let 
Uc  —  U(xc)  £  £(Rr,  Rn)  define  an  orthonormal  basis  of  ker  DG2(xc).  Then  the  implicit 
function  theorem  applied  to  the  equation 

(4.5  )  G2(xc  +  XJcy  +  DG2(xc)'^'  z)  =  0,  y  G  Rr,  z  £  Rp 

guarantees  the  existence  of  open  neighborhoods  Uc  of  the  origin  of  Rr  and  Vc  C  R"  of 
xc  such  that  for  any  y  £  Z7C  there  exists  exactly  one  solution  z  of  (4.5)  with  xc  +  Ucy  + 
DG2(xc)tz  £  Vc  and  that  the  mapping  ip  :  Uc  i— ►  Rp,  il'(y)  =  z  is  of  class  C1  on  Uc. 
Evidently,  we  have  V>(0)  =  0  and  Dip( 0)  =  0  and 

(4.6)  4>  :UC  Rn,  $(y)  =  xc  +  Ucy  +  DG2{xc)Tip(y),  VyeUc, 

is  a  diffeomorphism  from  Uc  onto  N  D  Vc-  In  other  words.  -1  is  a  chart  of  N  at  xc  and 
we  call  $  a  tangential  local  coordinate  map  at  xc. 

As  in  [RhSS],  by  shrinking  if  necessary  the  neighborhoods  Uc  and  Vc,  we  can  extend 
Uc  =  U{xc)  to  a  moving  frame  on  Vc;  that  is,  to  a  C^-mapping  U  :  Vc  — *  £(Rr,  R")  such 
that  the  columns  of  U(x)  form  an  orthonormal  basis  of  ker  DG 2(x)  for  each  x  £  Vc. 

Then,  for  y  €UC  and  x  =  $(y)  we  have  ( x,p )  £  M  exactljr  if 
B(y)p  =  H(y) 

(4.7)  B(y)  =  A\($(y))D$(y),  H(y)  =  G,(*(y)). 

Hence,  if  xc  £  W  then  necessarily  rank  Ai(x)  =  r  for  x  £  Rn  in  some  neighborhood  of  xc 
and  thus,  by  restricting  again,  if  needed,  the  neighborhoods  Uc  and  Vr  we  find  that  (4.7) 
represents  for  y  £  Uc  the  reduction  of  (4.1)  locally  near  xc. 

For  the  computation  we  need  to  be  able  to  evaluate  $(y)  and  D$(y)  for  y  £  Uc.  There 
are  various  possibilities  for  computing  x  —  $(y)  for  given  y.  For  example,  as  discussed  in 
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[RhSS],  we  may  use  the  QR-factorization 


R 

0 


(4.8)  DG2{xc)t  =  (QuQ2) 

where  rge  Q2  =  ker  DG2( xc)T  and  then  set  U°  =  Q2.  Now,  for  any  given  y  €  Rr  with 
sufficiently  small  norm  the  chord  Newton  algorithm 

Eval  $:  Input:  z,  R,  Qx 
while  ‘no  convergence’ 

solve  RT z  =  G2(x)  for  2; 
set  x  x  —  Q\Z\ 

Output:  x. 

converges  to  $(y)  €.  N  and  hence  implements  the  tangential  coordinate  system. 

For  the  computation  of  D$(y)  at  any  y  €  Uc  note  that 

(4.9)  ( Uc)TD$(y )  =  (UC)TUC  +  ( Uc)TDG2(xc)TDxP(y )  =  Ir. 

Moreover,  because  of  DG2($(y))D$(y)  =  0  it  follows  that  D$(y)  =  U($(y))K(y)  for  some 
nonsingular  I\(y)  €  £( Rr).  Together  with  (4.9)  this  implies  that  I\(y)  =  U($(y))TD$(y) 
and  therefore  that 

(4.10)  D*(y)  =  U(${y)){Uc)TU{${y)))~\ 

Clearly,  since  U  :  Vc  — »  £(3Rr,  Rn)  is  of  class  C1,  the  same  holds  for  £>#.  But  it  turns 
out  that  we  do  not  need  U(x)  to  be  a  C1 -moving  frame  on  a  neighborhood  of  xc  on  N.  In 
fact,  suppose  that  U(x)  6  £(Rr,  R^)  represents  for  z  =  ^(j/),  y  €  £/c,  an  arbitrarjr  basis 
matrix  of  ker  DG2(x )  and  that  U(xc)  =  U(xc).  Then  we  have  U(x )  =  U(x)Q(x )  with 
some  nonsingular  Q(x)  £  £(Rr)  and  hence 

D$(y)  =  U(x)((Uc)TU(x))~1  =  U(x)Q(x)((Uc)TU(x)Q(x))~1 

(4.11)  =  U(x)((Uc)TU(x))~1 . 
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While  for  points  x  in  a  neighborhood  of  a  fixed  point  xc  the  particular  choice  of  the 
basis  matrix  U(x )  does  not  matter,  the  orientation  of  U(xc)  =  Uc  does  play  a  role  when 
we  move  from  xc  to  another  “current”  point  xc.  The  compatibility  condition  for  charts  on 
a  C 1  manifold  requires  at  least  that  Uc  =  U(xc )  tends  to  Uc  =  U(xc)  when  xc  converges 
to  xc.  This  can  be  guaranteed  by  applying  the  moving  frame  algorithm  of  [RhSS]  in  the 
construction  of  the  new  basis  Uc.  However,  when  (4.1)  is  a  nonsingular  DAE  of  index  1 
and  hence  when  one  of  the  augmentation  procedures  of  section  2  can  be  applied  to  the 
reduced  system  (4.7),  then  it  suffices  to  ensure  that  both  bases  have  the  same  orientation; 
that  is,  that 

(4.12)  det  (VC)TUC  >  0. 


This  is  certainly  ensured  by  the  moving  frame  algorithm.  But  in  practice,  it  was  found 
advantageous  to  apply  a  simpler  heuristic  procedure.  Let  uc  =  £7C(1, 1, . . .  ,1)T  be  the 
“diagonal”  vector  of  the  positive  octant  of  the  basis  Uc.  Suppose  that  a  new  basis  matrix 
Uc  has  been  computed  with  the  columns  u1, . . .  ,ur.  Then  we  replace  the  vector  u'  by  —  ul 
whenever  ( uc)T u 1  <  0  for  any  i  =  1, . . .  ,  r.  This  certainly  guarantees  that  (4.12)  holds  for 
the  modified  basis. 

As  indicated,  if  (4.1)  is  a  nonsingular  DAE  of  index  1  then  one  of  the  augmentation 
procedures  of  section  2  will  be  applied  to  the  reduced  ODE  (4.7).  In  practice,  it  is  useful 
to  work  with  a  larger  augmentation  (2.14)  rather  than  with  (2.8)  in  order  for  the  process 
to  function  also  near  higher  order  singularities  than  just  standard  impasse  points.  For 
this  the  augmentation  (2.22)  with  q  —  2  was  chosen  for  which  v(y )  and  7 (y)  are  easily 
determined  by  (2.23). 

For  the  solution  of  (4.1)  subject  to  some  initial  condition  x(0)  =  x°  €  N  the  augmenta¬ 
tion  is  constructed  at  certain  computed  points  xc  along  the  trajectory.  At  these  points  the 
local  coordinate  map  (4.6)  is  obtained  and  hence  the  reduced  system  (4.7)  has  the  form 

Bcp  =  Hc ,  Bc  =  A,(xc)Uc,  Hc  =  Gn(xc) 
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There  are  various  ways  for  computing  suitable  vectors  b,  ci,c2  to  ensure  that  the  matrix 
of  the  augmented  system  (2.22);  that  is,  here 


Bc 

-Hc 

e\ 

(v' 

v2 

T 

C1 

0 

0 

Wi 

1 

W2 

rT 

C2 

0 

0/ 

L 
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(4.13) 


is  nonsingular.  For  relatively  small  dimensions  r  we  may  use,  for  instance,  the  singular 
value  decomposition 


(4.14) 


V?(B',-H')VR  =  (S,  0), 

VL  €  Rrxr,  Fs€Kr+1><r+1,  E  =  diag(<7, . a,). 


Let  e™,  k  =  1, . . .  ,  m,  denote  the  natural  unit-basis  vectors  of  Rm.  Then,  with 

=  VRe^  =  (X'2\ 

\^1  /  \^2/ 

we  choose  the  augmenting  vectors 

U-15)  e  =  ±VLe;,cl  =  -!!\,c1  =  r!^r, 

llUl|l2  IIU2||2 

where  the  sign  of  e  will  be  addressed  shortly.  Under  the  assumption  that  rank  (J5C,  -Hc)  > 
r  -  1,  the  matrix  of  (4.13)  is  nonsingular  if  and  only  if  ^  1.  In  fact,  for  any  null- 

vector  w  =  {qT,t,r])T  €  Rr+2  of  (4.13)  we  have 


(4.16) 


<7,5,  =  0,  *  =  !,...  ,  r  —  1,  (7rqr  -b  77  =  0,  for 


Hence,  in  the  case  of  rank  ( Bc,-Hc )  =  r  we  see  that  aT  >  0  and  therefore  that  77  =  0 
whence  q  =  Art},  £  =  Au>i  and  0  =  cjq  =  A;  that  is,  w  =  0.  When  the  rank  of  (Bc,  —Hc) 
is  r  —  1  then  (4.16)  implies  that  q  =  qrerr  and  therefore  that 

0  =  ujq  =  qrujii2  + 

0  =  u2q  =  qTu\u2  + 
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where,  because  of  the  orthonormality  of  the  vectors  Vjte^X  and  Vntrr+X  the  determinant 
of  the  matrix  equals  4-  —  1.  Hence,  if  this  determinant  is  nonzero  then  we  have 

again  w  =  0  while,  it  is  readily  seen,  that  in  the  case  of  a  zero  determinant  the  matrix 
may  indeed  be  singular.  Note  that  rank  ( Bc,—Hc )  >  r  —  1  obviously  holds  not  only  at 
standard  impasse  points  xc. 

For  the  choice  of  the  direction  of  the  vector  e  of  (4.15)  suppose  that  Bc  is  nonsingular; 
that  is,  that  we  are  not  exactly  at  a  singular  point  of  the  reduced  system.  Then  a  block 
LU-factorization  of  the  matrix  of  (4.13)  shows  that 

/  Bc  -Hc  e  \ 

(4.17)  det  cj  0  0  =  det  Bcdet  (cuc2)T(Bc)~1(Hc,-e). 

\  cj  0  0  / 

The  solution  of  the  augmented  system  (4.13)  gives 


(4.20) 


/ Bc  -Hc  e\ 

sign  det  cj  0  0  =  sign  det  Bc  sign  det  7' 

\cl  0  0/ 


In  line  with  the  theory  of  section  2  we  choose  the  augmenting  vector  e  such  that  the  left 
side  of  (4.20)  remains  positive.  In  other  words,  we  use  the  factorization  of  the  matrix  of  the 
augmented  system  (4.13)  to  monitor  the  sign  of  its  determinant  and  replace  the  computed 
zj  and  Z2  by  their  negative  values  when  that  determinant  is  negative. 


23 


For  the  computation  of  the  solution  of  (4.1)  with  the  initial  condition  .r(0)  =  .r°  we  use 
a  standard  explicit  Runge-Kutta  solver  such  as  RKF-45.  The  algorithm  for  one  step  along 
the  solution  of  (4.1)  then  has  the  general  form: 

SingDAE:  Input:  Current  point  xc,  tolerance  #o/, 

suggested  step  h;  minimal  step  hmin,  step  counter  k; 

1  Compute  the  QR-factorization  (4.S),  set  Uc  —  Q 2  and  yc  =  0; 

2  With  B{yc)  =  A\(xc)Uc  compute  the  augmentation  vectors  (4.13); 

3  Solve  the  augmented  system  (2.17)  to  obtain  ~,c  =  7(xc); 

4.  Take  a  Runge  Kutta  step  for  (2.12)  with  yc.  h.  to 1  to  obtain  yn ,  hntw : 

5.  If  ‘Step  not  accepted’  then 

replace  h  by  hj 2  and  for  |/t|  >  hmin  go  to  2,  otherwise  error  return: 

6.  Step  acceptable:  Call  Eval-#  to  obtain  xn  =  #(y")  G  N\ 

S.  Use  (4.10)  to  compute  £>#(yn)  and  B(yn)  =  Ai(xn)D$(yn); 

9.  Solve  the  augmented  system  (2.17)  to  obtain  7" 

10.  Output  k  =  k  +  1,  h  =  hnew,  x",  7C,  7"; 

11.  If ‘sign  7C  7^  sign  7" ’then 

singular  point  passed,  if  desired,  call  root  finder  to  compute  the  point. 

The  Runge  Kutta  solver  in  step  4.  requires  a  subroutine  for  computing  the  right  hand 
v(y)  of  the  reduced  explicit  equation  (4.7).  This  is  handled  by  a  subroutine  of  the  form: 

Eval  v:  Input :  Qi ,  and  R ,  y; 

1.  Use  Eval-#  to  compute  x  =  #(y); 

2.  Use  (4.10)  to  compute  D#(y)  and  B(y)  =  .4j (x)£>#(y); 

3.  Form  the  augmented  system  (2.17)  and  solve; 

4.  Use  (2.18)  to  compute  v(y)  and  7(y); 

5.  Output:  r(y).  7(y). 


24 


When  7  has  different  sign  at  two  consecutive  points,  say,  xk  and  xk+1 .  indicating  the 
presence  of  a  singular  point,  then  a  root  finder  can  be  applied  to  determine  the  step  h  from 
xk  which  gives  y(y)  =  0  and  hence  which  provides  the  singular  point.  For  this  a  simple 
algorithm  of  the  Dekker-Brent  type  (see  [B73])  has  been  used.  If  a  standard  impasse  point 
has  been  passed  then  xk  and  xk+1  lie  on  trajectories  with  opposite  orientation.  Hence,  in 
order  to  proceed  from  rfe+1,  we  have  to  change  the  sign  of  the  step  h.  At  higher  order 
singular  points  this  may  or  may  not  be  required. 

The  routine  ‘SingDae’  computes  a  new  augmentation  for  each  Runge  Ivutta  step.  This 
has  been  found  adequate  in  smaller  applications.  However,  for  larger  problems  it  is  desir¬ 
able  to  retain  the  same  augmentation  for  several  Runge  Kutta  steps.  For  this  the  condition 
of  the  linear  augmented  system  may  be  monitored.  This  is  easilv  done  in  step  4  of  ‘Eval 
v'  or  in  step  9.  of  ‘SinDae  Step".  An  additional  check  is  the  rate  of  convergence  of  the 
chord-Newton  process  of  ‘Eval  <J>\  Whenever  one  of  the?°  tests  indicates  that  the  augmen¬ 
tation  has  become  unreliable  then  the  output  of  ‘SingDae’  is  not  accepted  and  the  routine 
is  restarted  with  a  new  augmentation. 

There  should  be  no  need  to  enter  into  the  details  of  such  a  modification  of  the  process. 
The  approach  is  similar  to  that  employed  in  [PR91], 

5.  Numerical  Examples. 

As  noted  in  the  introduction,  impasse  points  for  DAE’s  arise  frequently  in  nonlinear 
circuit  problems.  As  an  example  we  consider  a  simple  circuit  consisting  of  a  nonlinear 
resistor,  linear  capacitor,  and  linear  inductor  in  parallel.  The  characteristic  of  the  resistor 
is  given  by  u  =  7  +  i2  where  i  and  u  denote  the  corresponding  branch  current  and  voltage 
drop,  respectively.  The  example  was  considered  earlier  by  F.  Tokens  (see  [T76])  and  again 
in  [RR92]  and  is  modelled  by  the  DAE 

0  0  10 
0  0  0  1 
0  0  0  0 
0  0  0  0 
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Here  Xj  =  ij,  j  —  1,2, 3,  are  the  currents  in  the  three  branches,  x±  =  u  is  the  voltage  drop, 
and,  for  simplicity,  the  capacitance  and  inductance  were  normalized  to  one.  It  is  readily 
verified  that  (5.1)  is  a  DAE  on  all  of  R4  with  r  =  2  in  the  sense  of  Definition  3.1. 

In  the  notation  of  (4.1)  the  constraint  manifold  (4.2)  is  here  given  by 


(5.2) 


Ar  =  {x  e  R4  :  G2(x)  =  0},  G2(x) 


/xj  +  X2  +  X3 

\  x4  -  7  -  x\ 


In  this  case  we  can  define  the  following  global  coordinate  mapping  on  N 


(5.3) 


:R2 


-  Ar,  *(i/)  = 


/  y'  \ 

y-i 

~(yi  +  Vi) 

V  7  +  y? 


and  hence  the  system 


(5.4) 


3  V  2z/!  0  /  \  y2 


cor.  itutes  a  (global)  reduction  of  (5.1). 

Clearly,  (5.4)  is  a  nonsingular  ODE  for  all  y  £  R2  with  yx  ^  0.  Moreover,  all  points 
of  {y  €  R2  :  yi  =  0,  y  ^  0}  are  standard  singular  points  of  (5.4)  while  at  y  =  0  the 
condition  (2.3)  is  violated.  Thus  (5.1)  has  the  one-dimensional  submanifold  {.r  6  R4;  x  = 
(0,£,  — ^ i  £  €  R,  i  7^  0}  of  standard  impasse  points  while  x*  =  (0,0,0, 7)r  is  a  higher 
singularity. 

A  closer  analysis  of  the  example  in  [RR92]  shows  that  the  singularity  at  x*  has  a 
different  character  for  the  four  cases  (i)  7  >  1/8,  (ii)  1/8  >  7  >  0,  (iii)  7  =  0,  (iv)  7  <  0. 

In  particular,  for  (ii)  -  (iv)  the  point  x*  is  a  funnel  point  (see  [T76]). 

Figures  1  -  3  show  computed  results  obtained  with  the  algorithm  of  Section  4  for  the 
three  values  7  =  1.0,  0.0625,-1.0.  In  each  case  the  impasse  points  x  =  (0.  £.  —  £,  7  )r 

for  £  >  0  are  accessible  while  those  for  £  <  0  are  inaccessible.  Moreover,  for  7  =  —1.0 
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the  funnel-nature  of  x*  is  clearly  visible,  but  recall  that  at  x*  the  derivative  x  becomes 
infinite  and  hence  no  C'-solution  passes  through  that  point.  The  results  certainly  show 
that  already  very  simple  circuits  may  have  a  relatively  complex  singularity  behavior. 

As  a  second  example  we  consider  a  two-phase  plug  flow  problem  described  by  Byrne 
and  Hindmarsh  ([BHS7],  see  also  [HLR89]).  The  equations  are  given  in  the  form: 


«)j£p(R  -  ^)2V^(2.5log[y^v/=p7  _  5]  +  10.5) 

Po 


R 


-bQcO--£QcO(l-b)  =  0, 

fpRy 


2n\l 5^<(2.W*  ~  l-25y2 ) log[ y  — -  5] 


+  3Ry  —  2.125y2  -  13 .6Rp 


pR 


)~Qa  =  0, 


where  P'  is  the  time-derivative  of  P.  With  xj  =  P,  x2  =  \/y/—  P\  and  .r3  =  y\J—P'  these 
equations  become 


x\x\  +1  =  0 

(5.5)  xi(4.2  +  log(cix3  -  5))(1  -  x2x3)2  -  c2x2(bx i  +  c4(l  -  b))  =  0 


a-3((2  -  x2x3)log(c1x3  -  5)  +  (2.4  -  1.7x2x3)  -  c3  =  0. 


Here,  after  scaling  the  time  by  a  factor  10  ,  the  constants  are 


R  pR 

Cl  =  7Vl-’  C2  = 


pQcO 


2.57rpf?ci  ’ 


C3  = 


10.8S  OApQg 

7T j.lRC] 


c4  —  10  Pq 


As  in  [HLR89]  we  used  the  values 


R  =  45.72,  p  =  0.814,  p  =  0.098,  b  =  0.345, 

Qco  =  1-7153  x  106,  Qa  =  3.027  x  105,  P0  =  1.378  x  108 


for  which 


ri  -  2012.47,  c2  =  19.7157,  c3  =  3.48464,  c4  =  13.7S 
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It  turns  out  that  the  system  (5.5)  has  a  singularity  at  x*  =  {0,0, x^)7,  with  = 
0.236849.  More  specifically,  the  system  has  index  one  in  the  sense  of  Definition  3.2  and  for 
the  reduced  ODE.  x *  is  a  basic  1-singular  point.  But,  since  the  form  (2.6)  is  zero  at  x* , 
this  point  is  not  a  standard  singular  point.  In  fact,  x *  corresponds  to  a  hysteresis  point 
for  stationary  equations. 

The  code  was  applied  to  the  system  (5.5)  starting  from  the  point 

y°  =(13.78,  0.42256012,  0.24921339)r 

given  in  [HLR89].  The  resulting  trajectory  if  shown  in  Figure  4.  The  singularity  y *  is 
reached  for  t*  =  1.0958  where  the  derivative  y'  becomes  infinite.  At  the  singular  point 
the  flow  is  choked  and  hence  the  part  of  the  trajectory  beyond  that  point  is  physically 
meaningless.  But  it  is  noteworthy  that  the  code  has  no  difficulty  in  passing  through  that 
point. 
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*3  -2-10  1 

gamma  =  -i.o 

Figure  3. 


t*  t 


Figure  4. 
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